Dispersive qubit readout using a thermal cavity

Ondrej Cernotik

23 February 2015

The evolution of a qubit dispersively coupled to a thermal cavity that is countinuously monitored is given by the stochastic master equation $$ \mathrm{d}\rho = -i\chi[\sigma_z r_\phi, \rho]\mathrm{d}t+\kappa\{(\bar{n}+1)\mathcal{D}[a]+\bar{n}\mathcal{D}[a^\dagger]\}\rho\mathrm{d}t+\sqrt{\frac{\kappa}{2\bar{n}+1}}\mathcal{H}[(\bar{n}+1)a-\bar{n}a^\dagger]\rho\mathrm{d}W, $$ where $\chi$ is the qubit-cavity coupling, $r_\phi = (ae^{i\phi}+a^\dagger e^{-i\phi})/\sqrt{2}$ is a quadrature of the cavity field which couples to an environment with $\bar{n}$ thermal photons at a rate $\kappa$, $\mathcal{D}[O]\rho = O\rho O^\dagger-\frac{1}{2}(O^\dagger O\rho+\rho O^\dagger O)$ is a Lindblad term, and $\mathcal{H}[O]\rho = (O-\langle O\rangle)\rho+\mathrm{H.c.}$ is the measurment term. Since the QuTiP stochastic solver [1,2] treats the measurement operators also as Lindblad operators, i.e., corresponding to decay, it is necessary to rewrite the Lindblad terms in the stochastic master equation using the measurement operator. We have $$ \kappa(\bar{n}+1)\mathcal{D}[a]\rho+\kappa\bar{n}\mathcal{D}[a^\dagger]\rho = \frac{\kappa}{2\bar{n}+1}\mathcal{D}[(bar{n}+1)a-\bar{n}a^\dagger]\rho +\frac{2\kappa\bar{n}(\bar{n}+1)}{2\bar{n}+1}\mathcal{D}[q]\rho $$ with $q = (a+a^\dagger)/\sqrt{2}$ being the amplitude quadrature of the cavity field.

When the cavity evolves on a time scale much faster than the qubit-cavity coupling, $\kappa\gg\chi$, its dynamics can be adiabatically eliminated from the equation of motion. Using Gaussian adiabatic elimination [3], the qubit dynamics is described by the equation $$ \mathrm{d}\rho_q = \frac{2\chi^2}{\kappa}(2\bar{n}+1)\mathcal{D}[\sigma_z]\rho_q\mathrm{d}t+\sqrt{\frac{2\chi^2}{\kappa(2\bar{n}+1)}}\mathcal{H}[-i(2\bar{n}\cos\phi+e^{-i\phi})\sigma_z]\rho_q\mathrm{d}W. $$ We compare this method of adiabatic elimination with results obtained by density operator expansion [4] (this method is, in fact, used only to obtain the measurement term; the decay has been obtained using projection operator method [5] to account for thermal noise), which gives the equation $$ \mathrm{d}\rho_q = \frac{2\chi^2}{\kappa}(2\bar{n}+1)\mathcal{D}[\sigma_z]\rho_q\mathrm{d}t+\sqrt{\frac{2\chi^2}{\kappa}}\mathcal{H}[e^{-i(\phi+\pi/2)}\sigma_z]\rho_q\mathrm{d}W; $$ see Ref. [3] for details.

The elimination methods and the exact solution are quantitatively compared using the trace distance of the two solutions. For a qubit, the trace distance of two density matrices is conveniently expressed using the expecation values of the Pauli operators $$ D(\rho_1, \rho_2) = \frac{1}{2}\sqrt{\sum_{i\in\{x,y,z\}}(\langle\sigma_i^1\rangle-\langle\sigma_i^2\rangle)^2}, $$ with $\langle\sigma_i^j\rangle = \mathrm{tr}(\rho_j\sigma_i)$. We use ensemble averaging to obtain an average trace distance for each elimination method as a function of time; additional time averaging gives a single figure of merit for the performance of the elimination methods.

References

[1] J. R. Johansson, P. D. Nation, and F. Nori, Computer Physics Communications 183, 1760 (2012).

[2] J. R. Johansson, P. D. Nation, and F. Nori, Computer Physics Communications 184, 1234 (2013).

[3] O. Cernotik, D. Vasilyev, and K. Hammerer, arXiv:1503.07484.

[4] A. Doherty and K. Jacobs, Physical Review A 60, 2700 (1999).

[5] C. Gardiner and P. Zoller, Quantum Noise: A Handbook of Markovian and Non-Markovian Quantum Stochastic Methods with Applications to Quantum Optics (Springer, 2004).


In [1]:
%matplotlib inline


Vendor:  Continuum Analytics, Inc.
Package: mkl
Message: trial mode expires in 29 days

In [2]:
from qutip import qeye, destroy, sigmax, sigmay, sigmaz, smesolve, ket2dm, basis, tensor, thermal_dm
from numpy import linspace, cos, sin, pi, sqrt, exp
from matplotlib.pyplot import subplots, plot

In [3]:
def trace_dist(chi, phi, kappa, nbar, Ntraj, NFock, Nstep, Nsub, tmax, rhoq0):
    #***** function trace_dist *****
    #This function calculates the average trace distance for Gaussian adiabatic
    #elimination and density operator expansion as functions of time.
    #Input parameters:
    #      chi: qubit-cavity coupling,
    #      phi: phase of the field quadrature in the interaction,
    #      kappa: intrinsic measurement (and decay) rate,
    #      nbar: thermal occupation,
    #      Ntraj: number of generated quantum trajectories,
    #      NFock: Fock number cutoff,
    #      Nstep: number of steps for evolution,
    #      Nsub: number of substeps for the stochastic solver,
    #      tmax: maximum time,
    #      rhoq0: initial qubit state.
    #Returns:
    #      tlist: list of times,
    #      D_gae: trace distance for Gaussian adiabatic elimination,
    #      D_doe: trace distance for density operator expansion.
    
    #***** System operators *****
    tlist = linspace(0,tmax,Nstep)
    Id2 = qeye(2)
    IdN = qeye(NFock) # identity operators on the qubit and the cavit field
    a = tensor(Id2,destroy(NFock)) # field annihilation operator
    q = (a+a.dag())/sqrt(2.)
    p = -1j*(a-a.dag())/sqrt(2.) # field quadratures
    sx = tensor(sigmax(),IdN)
    sy = tensor(sigmay(),IdN)
    sz = tensor(sigmaz(),IdN) # Pauli matrices
    
    #***** Full dynamics *****
    H_full = chi*sz*(q*cos(phi)-p*sin(phi)) # Hamiltonian
    c_op_full = [sqrt(2*kappa*nbar*(nbar+1)/(2*nbar+1))*q] # jump operators
    sc_op_full = [sqrt(kappa/(2*nbar+1))*((nbar+1)*a-nbar*a.dag())] # measurement and jump operators
    e_op_full = [sx, sy, sz] # operators whose evolution we are interested in
    rho0 = tensor(rhoq0, thermal_dm(NFock, nbar)) # initial state
    
    #***** Gaussian adiabatic elimination *****
    H_gae = 0*Id2
    c_op_gae = [sqrt(2*chi*chi*(2*nbar+1)/kappa-(2*chi*chi/(kappa*(2*nbar+1)))*(2*nbar*cos(phi)+exp(-1j*phi))*(2*nbar*cos(phi)+exp(1j*phi)))*sigmaz()]
    sc_op_gae = [-1j*sqrt(2*chi*chi/(kappa*(2*nbar+1)))*(2*nbar*cos(phi)+exp(-1j*phi))*sigmaz()]
    e_op_gae = [sigmax(), sigmay(), sigmaz()]
    
    #***** Density operator expansion *****
    H_doe = H_gae
    c_op_doe = [sqrt(4*chi*chi*nbar/kappa)*sigmaz()]
    sc_op_doe = [sqrt(2*chi*chi/kappa)*exp(-1j*(phi+pi/2.))*sigmaz()]
    e_op_doe = e_op_gae

    #***** Generating trajectories *****
    D_gae = 0
    D_doe = 0
    
    for i in xrange(0,Ntraj):
        print i
        sol_full = smesolve(H_full, rho0, tlist, c_op_full, sc_op_full, e_op_full, nsubsteps=Nsub, method='homodyne', 
                            solver='fast-milstein', store_measurement=True)
        sol_gae = smesolve(H_gae, rhoq0, tlist, c_op_gae, sc_op_gae, e_op_gae, nsubsteps=Nsub, method='homodyne', 
                           solver='fast-milstein', noise=sol_full.noise)
        sol_doe = smesolve(H_doe, rhoq0, tlist, c_op_doe, sc_op_doe, e_op_doe, nsubsteps=Nsub, method='homodyne', 
                           solver='fast-milstein', noise=sol_full.noise)
        s1x = sol_full.expect[0]
        s1y = sol_full.expect[1]
        s1z = sol_full.expect[2]
        s2x = sol_gae.expect[0]
        s2y = sol_gae.expect[1]
        s2z = sol_gae.expect[2]
        s3x = sol_doe.expect[0]
        s3y = sol_doe.expect[1]
        s3z = sol_doe.expect[2]
        
        D_gae = D_gae+sqrt((s1x-s2x)**2+(s1y-s2y)**2+(s1z-s2z)**2)/2.
        D_doe = D_doe+sqrt((s1x-s3x)**2+(s1y-s3y)**2+(s1z-s3z)**2)/2.

    D_gae = D_gae/Ntraj
    D_doe = D_doe/Ntraj
    
    return tlist, D_gae, D_doe

Calculating average trace distance -- ensemble averaging


In [4]:
chi = 0.1
phi = pi/2.
kappa = 1.
nbar = 0.5
Ntraj = 20
NFock = 8
Nsteps = 1000
Nsub = 250
tmax = 50

t, D_gae, D_doe = trace_dist(chi, phi, kappa, nbar, Ntraj, NFock, Nsteps, Nsub, tmax, ket2dm((basis(2,0)+basis(2,1)).unit()))


0
Total run time:  15.47s
Total run time:  10.19s
Total run time:   9.97s
1
Total run time:  15.13s
Total run time:   9.82s
Total run time:  10.01s
2
Total run time:  14.78s
Total run time:   9.80s
Total run time:   9.85s
3
Total run time:  14.93s
Total run time:   9.86s
Total run time:   9.98s
4
Total run time:  15.22s
Total run time:   9.82s
Total run time:   9.83s
5
Total run time:  15.06s
Total run time:   9.59s
Total run time:   9.49s
6
Total run time:  14.58s
Total run time:   9.71s
Total run time:   9.55s
7
Total run time:  14.51s
Total run time:   9.47s
Total run time:   9.74s
8
Total run time:  14.95s
Total run time:   9.94s
Total run time:   9.92s
9
Total run time:  15.00s
Total run time:   9.67s
Total run time:   9.84s
10
Total run time:  14.99s
Total run time:   9.55s
Total run time:   9.49s
11
Total run time:  14.63s
Total run time:   9.55s
Total run time:   9.99s
12
Total run time:  14.56s
Total run time:   9.53s
Total run time:   9.46s
13
Total run time:  14.63s
Total run time:   9.59s
Total run time:   9.53s
14
Total run time:  14.60s
Total run time:   9.79s
Total run time:   9.65s
15
Total run time:  14.91s
Total run time:   9.74s
Total run time:   9.60s
16
Total run time:  14.93s
Total run time:   9.84s
Total run time:   9.83s
17
Total run time:  15.72s
Total run time:   9.57s
Total run time:  10.08s
18
Total run time:  14.89s
Total run time:  10.67s
Total run time:  10.64s
19
Total run time:  15.26s
Total run time:  10.79s
Total run time:   9.87s

In [5]:
fig, axes = subplots(figsize=(12,8))

axes.plot(t, D_gae, label='Gaussian ad. elimination')
axes.plot(t, D_doe, label='Density op. expansion')
axes.set_xlabel('Time')
axes.set_ylabel('Trace distance')
axes.legend(loc=0)


Out[5]:
<matplotlib.legend.Legend at 0x109fceb10>

Calculating average trace distance -- time averaging

Erorr with Gaussian adiabatic elimination

In [6]:
D_gae_time = D_gae.sum()/Nsteps
print D_gae_time


0.0574992960051
Error with density operator expansion

In [7]:
D_doe_time = D_doe.sum()/Nsteps
print D_doe_time


0.130030346148